/*
 * CPU Reference Calculation for GPU Correctness Verification
 * Provides ground truth calculations for validating GPU kernel outputs
 */

#include <cmath>
#include <cfloat.h>
#include <cstdio.h>
#include <cstdlib.h>
#include <cstring.h>

// CPU reference implementation for vector addition
class CPUReferenceCalculator {
private:
    int array_size;
    int num_sub_tasks;
    float tolerance;

public:
    CPUReferenceCalculator(int array_size, int num_sub_tasks, float tolerance = 1e-6f)
        : array_size(array_size), num_sub_tasks(num_sub_tasks), tolerance(tolerance) {}

    // Simple vector addition (used for validation)
    void vectorAddCPU(const float* a, const float* b, float* c, int n) {
        for (int i = 0; i < n; i++) {
            c[i] = a[i] + b[i];
        }
    }

    // Sequential execution mode reference
    void sequentialExecution(const float** a_arrays, const float** b_arrays, 
                           float** c_arrays, int array_size, int num_sub_tasks) {
        for (int task = 0; task < num_sub_tasks; task++) {
            vectorAddCPU(a_arrays[task], b_arrays[task], c_arrays[task], array_size);
        }
    }

    // Masquerading execution mode reference
    void masqueradingExecution(const float** a_arrays, const float** b_arrays,
                             float** c_arrays, int array_size, int num_sub_tasks) {
        int total_elements = array_size * num_sub_tasks;
        
        for (int global_idx = 0; global_idx < total_elements; global_idx++) {
            int sub_task_id = global_idx / array_size;
            int local_idx = global_idx % array_size;
            
            if (sub_task_id < num_sub_tasks && local_idx < array_size) {
                c_arrays[sub_task_id][local_idx] = 
                    a_arrays[sub_task_id][local_idx] + b_arrays[sub_task_id][local_idx];
            }
        }
    }

    // Sub-task sequential execution reference
    void subTaskSequentialExecution(const float** a_arrays, const float** b_arrays,
                                  float** c_arrays, int array_size, int num_sub_tasks) {
        for (int task = 0; task < num_sub_tasks; task++) {
            vectorAddCPU(a_arrays[task], b_arrays[task], c_arrays[task], array_size);
        }
    }

    // Memory layout analysis functions
    void generateLinearPattern(float* array, int size, float start = 0.0f, float step = 1.0f) {
        for (int i = 0; i < size; i++) {
            array[i] = start + i * step;
        }
    }

    void generateInterleavedPattern(float* array, int size, int stride, int offset = 0) {
        for (int i = 0; i < size; i++) {
            array[i] = (float)(i * stride + offset);
        }
    }

    void generateBlockedPattern(float* array, int size, int block_size, int offset = 0) {
        for (int i = 0; i < size; i++) {
            int block = i / block_size;
            int pos_in_block = i % block_size;
            array[i] = (float)(block * block_size + pos_in_block + offset);
        }
    }

    // Edge case generation functions
    void generateZeroArrays(float* a, float* b, int size) {
        memset(a, 0, size * sizeof(float));
        memset(b, 0, size * sizeof(float));
    }

    void generateNaNArrays(float* a, float* b, int size) {
        for (int i = 0; i < size; i++) {
            if (i % 100 == 0) {  // Every 100th element
                a[i] = NAN;
                b[i] = 1.0f;
            } else {
                a[i] = (float)i;
                b[i] = (float)i;
            }
        }
    }

    void generateInfinityArrays(float* a, float* b, int size) {
        for (int i = 0; i < size; i++) {
            if (i % 50 == 0) {  // Every 50th element
                a[i] = INFINITY;
                b[i] = 1.0f;
            } else {
                a[i] = (float)i;
                b[i] = (float)i;
            }
        }
    }

    void generateDenormalArrays(float* a, float* b, int size) {
        for (int i = 0; i < size; i++) {
            if (i % 200 == 0) {  // Every 200th element
                a[i] = 1e-39f;  // Denormal value
                b[i] = 1e-39f;
            } else {
                a[i] = (float)i;
                b[i] = (float)i;
            }
        }
    }

    // Validation functions
    bool validateResults(const float* computed, const float* expected, int size, 
                        float tolerance = 1e-6f, int* error_count = nullptr) {
        int errors = 0;
        
        for (int i = 0; i < size; i++) {
            float diff = fabsf(computed[i] - expected[i]);
            
            // Handle special cases
            if (isnan(computed[i]) && isnan(expected[i])) {
                continue;  // Both NaN
            } else if (isinf(computed[i]) && isinf(expected[i])) {
                continue;  // Both infinity
            } else if (isnan(computed[i]) || isnan(expected[i])) {
                errors++;
                if (error_count) *error_count = errors;
                continue;
            } else if (isinf(computed[i]) || isinf(expected[i])) {
                errors++;
                if (error_count) *error_count = errors;
                continue;
            } else if (diff > tolerance) {
                errors++;
                if (error_count) *error_count = errors;
            }
        }
        
        if (error_count) *error_count = errors;
        return errors == 0;
    }

    // Detailed validation with error reporting
    bool validateWithReporting(const float* computed, const float* expected, int size,
                             float tolerance = 1e-6f, bool verbose = false) {
        int errors = 0;
        float max_error = 0.0f;
        int max_error_index = -1;
        
        for (int i = 0; i < size; i++) {
            float diff = fabsf(computed[i] - expected[i]);
            
            if (diff > tolerance || 
                (isnan(computed[i]) != isnan(expected[i])) ||
                (isinf(computed[i]) != isinf(expected[i]))) {
                
                errors++;
                if (diff > max_error) {
                    max_error = diff;
                    max_error_index = i;
                }
                
                if (verbose && errors <= 10) {  // Limit verbose output
                    printf("Error at index %d: computed=%.6f, expected=%.6f, diff=%.6f\n",
                           i, computed[i], expected[i], diff);
                }
            }
        }
        
        if (verbose) {
            printf("Validation: %d/%d elements correct, max error=%.6f at index %d\n",
                   size - errors, size, max_error, max_error_index);
        }
        
        return errors == 0;
    }

    // Memory layout verification
    void verifyMemoryLayout(const float** arrays, int num_arrays, int array_size,
                          const char* layout_name) {
        printf("Verifying memory layout: %s\n", layout_name);
        
        for (int i = 0; i < num_arrays; i++) {
            printf("Array %d: first=%.1f, last=%.1f\n", i, arrays[i][0], arrays[i][array_size-1]);
        }
    }

    // Performance model for expected results
    void calculateExpectedPerformanceModel(int array_size, int num_sub_tasks,
                                         float* expected_time, float* memory_bandwidth) {
        // Simple model: O(n) computation + memory transfer
        float computation_factor = 1e-9f;  // Rough estimate
        float memory_factor = 1e-9f;       // Memory bandwidth factor
        
        *expected_time = (array_size * num_sub_tasks * computation_factor) + 
                        (array_size * num_sub_tasks * 3 * memory_factor);
        *memory_bandwidth = (array_size * num_sub_tasks * 3 * sizeof(float)) / (*expected_time);
    }
};

// C-style interface functions
extern "C" {
    // Create reference calculator
    CPUReferenceCalculator* create_reference_calculator(int array_size, int num_sub_tasks, float tolerance) {
        return new CPUReferenceCalculator(array_size, num_sub_tasks, tolerance);
    }

    void destroy_reference_calculator(CPUReferenceCalculator* calculator) {
        delete calculator;
    }

    // CPU reference implementations
    void cpu_vector_add(const float* a, const float* b, float* c, int n) {
        for (int i = 0; i < n; i++) {
            c[i] = a[i] + b[i];
        }
    }

    void cpu_sequential_execution(const float** a_arrays, const float** b_arrays,
                                float** c_arrays, int array_size, int num_sub_tasks) {
        CPUReferenceCalculator calc(array_size, num_sub_tasks);
        calc.sequentialExecution(a_arrays, b_arrays, c_arrays, array_size, num_sub_tasks);
    }

    void cpu_masquerading_execution(const float** a_arrays, const float** b_arrays,
                                  float** c_arrays, int array_size, int num_sub_tasks) {
        CPUReferenceCalculator calc(array_size, num_sub_tasks);
        calc.masqueradingExecution(a_arrays, b_arrays, c_arrays, array_size, num_sub_tasks);
    }

    // Validation functions
    int validate_float_arrays(const float* computed, const float* expected, int size,
                            float tolerance, int* error_count) {
        CPUReferenceCalculator calc(size, 1, tolerance);
        return calc.validateResults(computed, expected, size, tolerance, error_count) ? 1 : 0;
    }

    void generate_test_data(float* a, float* b, int size, int pattern_type) {
        CPUReferenceCalculator calc(size, 1);
        
        switch (pattern_type) {
            case 0: // Linear
                calc.generateLinearPattern(a, size);
                calc.generateLinearPattern(b, size, 1.0f, 1.0f);
                break;
            case 1: // Interleaved
                calc.generateInterleavedPattern(a, size, 2);
                calc.generateInterleavedPattern(b, size, 2, 1);
                break;
            case 2: // Blocked
                calc.generateBlockedPattern(a, size, 1024);
                calc.generateBlockedPattern(b, size, 1024, 1);
                break;
            case 3: // Zeros
                calc.generateZeroArrays(a, b, size);
                break;
            case 4: // NaN test
                calc.generateNaNArrays(a, b, size);
                break;
            case 5: // Infinity test
                calc.generateInfinityArrays(a, b, size);
                break;
            default:
                calc.generateLinearPattern(a, size);
                calc.generateLinearPattern(b, size, 1.0f, 1.0f);
                break;
        }
    }
}